#pragma once
#include <iostream>
#include "GridLayer.h"
#include <GL/glew.h>
#include "CgProgram.h"
#include "RenderTexture.h"
#include "PenroseTiling.h"
#include "ParallelPoison.h"
#include "WangTiles.h"
#include "vmath.h"


// Class that implements fluid simulation on a 3D grid.
class FluidGridSimulator
{
public:

	FluidGridSimulator(int m, int n, int l, float sizeX = 1, float sizeY = 1, float sizeZ = 1, float posX = 0, float posY = 0, float posZ = 0);
	~FluidGridSimulator(void);

	//computes one simulation step of the fluid
	void update();
	//renders the fluid
	void render(int w, int h);
	//renders the bounding box
	void renderBoundingBox(float hx, float hy, float hz);
	//the fluid can be visualized in different ways, this function switches between them
	void switchDisplayMode()
	{
		displaymode = (displaymode + 1)%5;
		switch(displaymode)
		{
		case 0:
			sprintf(displayModeString,"NPR display");
			break;
		case 1:
			sprintf(displayModeString,"Real smoke display");
		case 2:
			sprintf(displayModeString,"Detected edges");
			break;
		case 3:
			sprintf(displayModeString,"Velocity Divergence");
			break;
		case 4:
			sprintf(displayModeString,"Pressure");
			break;
		case 5:
			sprintf(displayModeString,"Density");
			break;
		case 6:
			sprintf(displayModeString,"Vorticity");
			break;
		case 7:
			sprintf(displayModeString,"Vorticity Force");
			break;
		case 8:
			sprintf(displayModeString,"Temperature");
			break;
		case 9:
			sprintf(displayModeString, "Grid Point");
			break;
		}
	}
	int pointdisalg, edgemode, joinmode, gausssig, hatchmode, intensitymode;
	void switchPointDistributeAlg()
	{
		pointdisalg = (pointdisalg + 1)%3;
		genSamplingList();
	}
	void switchEdgeDetectionMode()
	{
		edgemode = (edgemode + 1)%2;
	}
	void switchHatchMode()
	{
		hatchmode = (hatchmode + 1)%2;
	}
	void switchIntensityMode()
	{
		intensitymode = (intensitymode + 1)%2;
	}
	void switchJoinEdgeMode()
	{
		joinmode = (joinmode + 1)%3;
	}
	//the npr smoke could be rendered with different kind of stroke
	void switchStrokeTexture()
	{
		strokeTexType = (strokeTexType + 1)%3;
		createColorBand();
	}
	void switchNPRStyle()
	{
		nprStyle = (nprStyle + 1)%3;
	}
	//pause or play smoke
	void setPause(bool value){pause = value;}
	bool getPause(){return pause;}
	//camera change?
	void setCameraChanged(bool value)
	{
		isCamChanged = value;
	}
	bool isCameraChanged(){return isCamChanged;}
	//turn on and off render based on density gradient
	void setDensityGradient(bool value){computeDensGrad = value;}
	bool getDensityGradient(){return computeDensGrad;}
	//turns on and off buoyancy force
	void setBuoyancyForces(bool value){computeBuoyancy = value;}
	bool getBuoyancyForces(){return computeBuoyancy;}
	// returns the display mode as a string
	char* getDisplayMode()
	{	
		return displayModeString;
	}
	//set fluid temperature to zero
	void resetTemperature();
	// set fluid density to zero
	void resetDensity();
	// set fluid velocity to zero
	void resetVelocity();
	//turns on and off early Z culling
	void setEarlyZCull(bool value){earlyZCull = value;}
	bool getEarlyZCull(){return earlyZCull;}
	// sets iteration count used in pressure computation
	int getIterationCount(){return iterationCount;}
	void setIterationCount(int count){if(count < 0) count = 0; iterationCount = count;}
	// sets fluid viscosity
	void setViscosity(float value){viscosity = value; if(viscosity < 0.001) viscosity = 0;}
	float getViscosity(){return viscosity;}
	//turn on/off show real smoke
	void setRealSmoke(bool value){realSmoke = value;}
	bool getRealSmoke(){return realSmoke;}
	// turns on/off vorticity confinement
	void setVorticityConfinement(bool value){computeVorticity = value;}
	bool getVorticityConfinement(){return computeVorticity;}
	// turns on/off bounding cube rendering
	void setRenderBoundingCube(bool value){renderBoundingCube = value;}
	bool getRenderBoundingCube(){return renderBoundingCube;}
	// sets slice count used at rendering
	void setRenderSliceCount(int count){renderSliceCount = count; if(renderSliceCount < 0) renderSliceCount = 0;}
	int getRenderSliceCount(){return renderSliceCount;}
	// tells the simulator to use the current external force texture in all frames
	void setLockExternalForces(bool value){lockExternalForces = value;}
	// add buoyancy force to the external forces
	void addBuoyancyForces(float dt);
	// add a force to the external forces
	void addExternalSpotSource(float* position = 0,
		float* direction = 0,
		float velocityRadius = 0,
		float densityRadius = 0,
		float tempRadius = 0,
		float velocityScale = 0,
		float densityScale = 0,
		float tempScale = 0);
	void addGaussianBall();
	void addExternalUniformSource(float density, float temp, float velocityX, float velocityY, float velocityZ);
	// returns the uppmost layer of the fluid
	void getUpperLayer(void* pixeldata);
	// returns the lowest layer of the fluid
	void getLowerLayer(void* pixeldata);
	// fills the upper bound layer with the given data
	void setUpperBoundLayer(void* pixeldata);
	// fills the lower bound layer with the given data
	void setLowerBoundLayer(void* pixeldata);
	// sets the type of simulation if distributed simulation is used
	void setType(int t){type = t;}
	// sets the flat3D texture that stores the voxelized obstacle geometry
	void setObstacleTexture(RenderTexture* rt)
	{
		obstacleTexture = rt;
		addObstacleInfluence = (rt != 0);
		if(!addObstacleInfluence)
		{
			glColorMask(false,false,false,true);

			pressure->startCopy();
			glClearColor(0,0,0,0);
			glClear(GL_COLOR_BUFFER_BIT);
			pressure->endCopy();

			glColorMask(true,true,true,true);			
		}
	}
	// set external forces to zero
	void clearExternalForces();

	RenderTexture* getExternalVelocityTexture(){return externalVelocity;}
	RenderTexture* getExternalDensityTexture(){return externalDensity;}
	GLuint getDensityTextureID(){return density->getTextureID();}
	void setVelocityBoundScale(float scale){velocityBoundScale = scale;}
	void setPressureBoundScale(float scale){pressureBoundScale = scale;}

	void setDisThres(float value){m_disthres = value;}
	float getDisThres(){ return m_disthres;}

protected:
	// slice count used at rendering
	int renderSliceCount;
	// the name of the current display mode
	char displayModeString[100];
	//turns on/off show realistic smoke
	bool realSmoke;
	// turns on/off vorticity confinement
	bool computeVorticity;
	// turns on/off bounding cube rendering during 3d display mode 
	bool renderBoundingCube;
	// iteration count used in pressure computation
	int iterationCount;
	//turns on and off early Z culling
	// voxels with low pressure values need fewer iteration steps than voxels with high pressure values
	// early Z cull uses depth test to limit the number of iterations for each voxel according to its pressure
	bool earlyZCull;
	//the fluid can be visualized in different ways
	int displaymode;
	// the type of simulation if distributed simulation is used
	// 0: no distribution
	// 1: the fluid grid-part simulated by this instance is at the bottom of the original grid (according to the z axis)
	// 2: the fluid grid-part simulated by this instance is somewhere in the middle of the original grid (according to the z axis)
	// 3: the fluid grid-part simulated by this instance is the topmost part of the original grid (according to the z axis)
	int type;
	// helper texture to speed up flat+D texture indexing
	GLuint layerOffsetTexture;
	// color texture for gradient rendering
	GLuint strokeTexture, nprcolorBand, halftoneBand, splatTexture, opaquediskTexture, contourstrokeTexture,
		colorstrokeBand, randTex;
	// external forces should be added during the simulation step
	bool externalForces;
	// buoyancy force should be added during the simulation step
	bool computeBuoyancy;
	//density gradient should be computed
	bool computeDensGrad;
	//pause or play
	bool pause;
	//camera changed?
	bool isCamChanged;
	// the current value of external forces should be used in all simulatin steps
	bool lockExternalForces;
	// only for testing, now it is always set to true
	bool clearPressurePerFrame;
	// only for testing, now it is always set to false
	bool computePressureBoundPerIter;
	// the higher this value the more resistent is the fluid to forces (it gets thicker)
	float viscosity;
	// grid resolution along the y axis
	int res_n;
	// grid resolution along the x axis
	int res_m;
	// grid resolution along the z axis
	int res_l;
	// number of slice tiles in a row of the flat3D texture
	int numTilesX;
	// number of slice tiles in a column of the flat3D texture
	int numTilesY;
	// grid size in the 3D space
	float sizeX;	
	float sizeY;
	float sizeZ;
	// grid position in the 3D space
	float posX;
	float posY;
	float posZ;
	// does obstacle influence need to be computed
	bool addObstacleInfluence;
	// Objects that describe one slice (layer) of the 3D grid
	GridLayer** gridLayers;
	// Upper and lower bounding layers
	BoundaryLayer* boundLayers[2];
	//light viewport
	int lw, lh;
	//distance threshold
	float m_disthres;

	float velocityBoundScale;
	float pressureBoundScale;

	GLdouble preModelViewProj[16];
	GLfloat* randIntensity;

	/////////////////////////////////////////////////////
	// ray casting textures
	RenderTexture* intersection;
	// buoyancy textures
	RenderTexture* buoyancyForce;
	// external temperature
	RenderTexture* externalTemp;
	// MacCormack scheme
	RenderTexture* phi_n_1_hat;
	RenderTexture* phi_n_hat;
	// light buffer
	RenderTexture* lightbuffer;

	RenderTexture *footprint, *marktex, *initpointtex, *tmpadvecttex, *tmpadvecttex_2v,
		*ndotetex, *depthtex;
	IteratedRenderTexture* resulttex, *imgbuff;
	IteratedRenderTexture *edgetex;

	/////////////////////////////////////////////////////
	// ray casting Cg GPU program

	//fragment and vertex program used in intersection texture building
	CgProgram* s_IntersecfrontProgram;
	CgProgram* s_IntersecbackProgram;
	CgProgram* s_LocateIntersectProgram;

	//vertex program used in grid points rendering
	CgProgram* s_gridpointProgram;
	//MacCormack scheme
	CgProgram* s_advectMacCormack;
	
	//program for smoothing
	CgProgram* s_smoothProgram, *s_smooth2DProgram;
	//program for real rendering
	CgProgram* s_realRenderProgram, *s_renderTexProgram;
	
	CgProgram* s_alphatestProgram;
	CgProgram* s_extractSilhouetteProgram;
	CgProgram* s_renderResultTex;

	///////////////////////////////////////////////////
	// Flat3D render textures

	// velocity field
	IteratedRenderTexture* velocity;
	// density field
	IteratedRenderTexture* density;
	//temperature field
	IteratedRenderTexture* temperature;
	// pressure field
	IteratedRenderTexture* pressure;
	// divergence of the velociy field
	RenderTexture* velocityDiv;
	// density field from external sources
	RenderTexture* externalDensity;
	// velocity field from external sources
	RenderTexture* externalVelocity;
	// vorticity of the velocity field
	IteratedRenderTexture* vorticity;
	// velocity field created by vorticity
	RenderTexture* vorticityForce;
	// voxelized obstacle geometry
	RenderTexture* obstacleTexture;

	/////////////////////////////////////////////////////
	// Cg GPU programs
	//fragment and vertex program for hatching
	CgProgram* s_hatchingProgram, *s_intertexturesProgram, *s_advectSampPoint, *s_distPointProgram,  *s_distPointProgram2v,
		*s_markRegionProgram, *s_stippleProgram, *s_rotateSampPoint, *s_solidcontourProgram;
	CgProgram* s_edgedetectProgram, *s_texcontourProgram, *s_texcontourVProgram, *s_NPRcolorizeProgram, 
		*s_toonshadeProgram;
	CgProgram* s_billboardProgram;

	// fragment and vertex program for light accumulation
	CgProgram* s_lightProgram;
	//fragment program compute density gradient
	CgProgram* s_denGraProgram, *s_denHessianProgram;
	// fragment program computes the buoyancy forces
	CgProgram* s_BuoyancyProgram;	
	// this fragment program renders a constant colored fullscreen quad
	CgProgram* s_ColoredProgram;
	// this fragment program renders a textured fullscreen quad
	CgProgram* s_TexturedProgram;
	// fragment program used during advection of the velocity and density fields in internal voxels
	CgProgram* s_AdvectProgram;
	// fragment program used during advection of the velocity and density fields in boundary layer voxels
	CgProgram* s_BoundaryLayerProgram;
	// fragment program used during advection of the velocity and density fields boundary edge voxels
	CgProgram* s_BoundaryEdgeProgram;
	// this fragment program adds external velocity to the external sources
	CgProgram* s_AddForceProgram;
	// this fragment program adds external density to the external sources
	CgProgram* s_AddDensityProgram;
	// this fragment program computes velocity divergence
	CgProgram* s_DivergenceProgram;
	// this fragment program computes one Jacobi iteration step, used in pressure calculation and velocity diffusion
	CgProgram* s_JacobiIterProgram;
	// this fragment program does the projectio step
	CgProgram* s_ProjectProgram;
	// this fragment program computes the vorticity of the velocity field
	CgProgram* s_VorticityProgram;
	// this fragment program computes the velocity caused by vorticity
	CgProgram* s_VorticityForceProgram;
	// // this fragment program is used to write pressure values to the depth buffer (used in early Z culling)
	CgProgram* s_WritePressureAsDepthProgram;
	// vertex program used by the above fragment programs
	CgProgram* s_RenderQuadVertexProgram;
	// fragment and vertex program used in 3D display mode
	CgProgram* s_VolumeRenderSliceFlat3DProgram;


	// Functions of each of the simulation steps
	void Advect(float dt, IteratedRenderTexture* media, IteratedRenderTexture* velocity, float borderScale);
	void AddExternalForces();
	void AddObstacleInfluence();
	void ComputeVorticity(float dt);
	void DiffuseVelocity(float dt);
	void ComputePressure();
	void Project();
	void ComputeBoundary(IteratedRenderTexture* media, float borderScale, bool updateUpperBound = true, bool updateLowerBound = true);

	// render grid point
	void renderGridPoint(int w, int h);
	// render NPR smoke
	void renderFlat3DSliceQuad(int w, int h, IRenderTarget* v, int isSmoothed = false);
	// render image buff to the screen
	void renderRealSmoke(int w, int h, int w0 = 0, int h0 = 0);
	void renderTextures(int w, int h, IRenderTarget* rd);
	// render the fluid with the use of the flat3D density texture - ray casting
	void renderFlat3D(int w, int h);
	// renders a fullscreen quad
	void renderFullscreenQuad();
	// renders a fullscreen textured quad
	void renderTexturedQuad(IRenderTarget* colorTex, float trScale = 1, float trTrans = 0, CgProgram* program = 0);
	// create the helper texure used in flat3D texture indexing
	void createLayerOffsetTexture();		
	//create a color band texture used in gradient rendering
	void createColorBand();
	//cacalate transformation matrix of the light
	void calculateMatrix(float* light, GLdouble* lmvp, GLdouble* emvp, GLdouble* hwmv, bool &f2b);

	//create billboard texture
	void createPuffBillBoard();

	void ContourRendering();
	void NPRcoloring();
	void Hatching();
	void Stippling();
	void RorateSmoke();
	void SolidContour();
	inline void SurfaceDetect();

	float haltonsequence(int id, int base)
	{
		int I = id;
		float H = 0.f;
		float half = 1.f / base;
		float digit;
		while (I > 0)
		{
			digit = I % base;
			H = H + digit * half;
			I = (I*1.f - digit) / base;
			half = half / base;
		}
		return H;
	}

	//Light pos
	float lightPos[3];

	//test temporal coherence
	float dt;
	IteratedRenderTexture* renderedDensity;
	RenderTexture* renderedCurvature, *hessian_r1, *hessian_r2, *hessian_r3;
	RenderTexture* pointCollection, *pointCollection_2v, *advectCollection, *advectCollection_2v;
	RenderTexture* billBoard;

	//acceleration
	int dispListId, samplingListId;

public:
	int strokeTexType, nprStyle, innerStrokeType;
	float densityHaltonN;
	void Fading();
	void smoothingSmokeData();
	void genSamplingList();

public:
	//parameters
	float dsampleX, dsampleY, m_sampleRatio;
	int vpW, vpH;
	float m_size[3], m_pos[3];
	int nSmooth;
	int detectMode;
	float densityThres;
	float intensityThres;
	ParallelPoison ps;
	WangTile wang;

	bool is2Dsmoothed;
	float recsize;

	Vector3f m_impulse;
	float m_time;

	float m_windsource[3];
	float m_winddir[3];
	bool m_isWind;
	float m_contourtexsize;
	float m_contourthres;
	int m_sclstep;

public: 
	float getDt(){ return dt;}
	void setLightSource(float x, float y, float z)
	{
		lightPos[0] = x;
		lightPos[1] = y;
		lightPos[2] = z;
	}
	
	void initializeImpulse()
	{
		m_impulse.x = res_m / 2.f;
		m_impulse.y = res_n / 2.f;
		m_impulse.z = 1;

		m_windsource[0] = res_m / 2.f;
		m_windsource[1] = res_n / 2.f;
		m_windsource[2] = 1.f;

		m_isWind = false;
		m_winddir[0] = 0;
		m_winddir[1] = 0;
		m_winddir[2] = 0;
	}

	void setWindSource(float dx, float dy, float dz, int mod)
	{
		m_windsource[0] += dx;
		m_windsource[1] -= dy;
		m_windsource[2] -= dz;

		printf("windxyz: %f %f %f %d\n", m_windsource[0], m_windsource[1], m_windsource[2], m_isWind);

		if (mod == GLUT_ACTIVE_SHIFT)
		{
			float l = sqrt(dx*dx + dy*dy + dz*dz);
			if (l > 0.5f)
			{
				m_isWind = true;
				m_winddir[0] = dx / l;
				m_winddir[1] = dy / l;
				m_winddir[2] = dz / l;
			}
			else
				m_isWind = false;
		}
		else
		{
			m_isWind = false;
			m_winddir[0] = 0;
			m_winddir[1] = 0;
			m_winddir[2] = 0;
		}
	}
};
